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Abstract. A model is presented which calculates the gas temperature and chemistry in the surface layers of 
flaring circumstellar disks using a code developed for photon-dominated regions. Special attention is given to the 
influence of dust settling. It is found that the gas temperature exceeds the dust temperature by up to several 
hundreds of Kelvins in the part of the disk that is optically thin to ultraviolet radiation, indicating that the 
common assumption that T gaa = Td uat is not valid throughout the disk. In the optically thick part, gas and 
dust are strongly coupled and the gas temperature equals the dust temperature. Dust settling has little effect 
on the chemistry in the disk, but increases the amount of hot gas deeper in the disk. The effects of the higher 
gas temperature on several emission lines arising in the surface layer are examined. The higher gas temperatures 
increase the intensities of molecular and fine-structure lines by up to an order of magnitude, and can also have 
an important effect on the line shapes. 
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1. Introduction 

Circumstellar disks play a crucial role in both star and 
planet formation. After protostars have formed, part of 
the remnant material from the parent cloud core can con- 
tinue to accrete by means of viscous processes in the disk 
llShu et alJl!987h . Disks are also the sites of planet for- 
mation, either through coagulation and accretion of dust 
grains or through gravitat ional instabilities in the disk 
l|Lissaueilll993l iBossI 12000) . To obtain detailed informa- 
tion on these processes, both the dust and the gas com- 
ponent of the disks have to be studied. Over the last 
decades, there have been numerous models of the dust 
emission from disks, with most of the recent work fo- 
cussed on flaring structures in which the disk surface in- 
tercepts a significant part of the stellar radiation out to 
large distances and is heated to much higher tempera- 



tures than the midp] 


ane (e.g. iKenvon & HartmannI 


1987 


Chiang & Goldreich 


Il997t iD'Alessio et all Il998l 


1999 


Bell et al.ll997l 


Dullemond et al.12002"). These models are 



appropriate to the later stages of the disk evolution when 
the mass accretion rate has dropped and the remnant en- 
velope has dispersed, so that heating by stellar radiation 
rather than viscous energy release dominates. Such models 
have been shown to reproduce well the observed spectral 
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energy distributions at mid- and far-infrared wavelengths 
for a large variety of disks around low- and intermediate- 
mass pre-main sequence stars. 

Observations of the gas in disks started with millime- 
ter interferometry data of the lowest transitions of the 
CO molecule (e.g. lKoerner fc SargentHl 995HDutrev et all 



Il996t iMannings fc Sargent] 1 199 71: iDartois et all 1200*^ 
but now also include submillimeter single-dish (e.g 
[ Kastner et al1ll997t iThi et all EoOll Ivan Zadelhoff et al 
| 200l|) and infrared (e.g lNaiita et alJl2003l iBrittain et al 
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2003J) data on higher excitation CO lines. Evidence for 
the presence of warm gas in disks also comes from 
near- and mid-infr ared and ultraviolet observations of the 

molecule (e .g. iHerczeg et all 120021: iBarv et al.l l2003t 
IThi et al"ll200l[) . Although emission from the hottest gas 
probed at near-infrared and ultraviolet wavelengths is 
thought to come primarily from a region within a few AU 
of the young stars, the longer wavelength data trace gas at 
larger distances from the star, >50 AU. Molecules other 
than CO are now also detected at (sub)millimcter wave- 
lengths, including CN, HCN, HCO+, CS and H?CO (e. 
i Dutrev et aD Il997t iKastner et ail Il997t iQi et all 1200 
IThi et alJl2004h . 

The emission from all of these gas tracers is deter- 
mined both by their chemistry and excitation, where 
the latter depends on the temperature and density 
structure of the disk. The chemistry in flaring disks 
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Fig. 1. The 1+1-D model. The disk is divided into annuli, 
each of which is treated as a 1-dimensional PDR problem 
for the chemistry and the cooling rates. 

has been studie d intensely in recent years by vari- 
ous groups (e.g. lAikawa et al. I Il997t IWillacv fc Langel 
l200a iMarkwick et all 120021: Ivan Zadelhoff et alJ l2003|) . 

whereas the density structure is constrained from verti- 
cal hydrostatic equilibrium model s of the dust disk as- 
suming a contant gas/ dust ratio (D'Aless io et alJ ll999: 
iDullemond et al1l2002l) . In all of these models, the gas 
temperature is assumed to be equal to the dust temper- 
ature. Th at this assumption doe s not always hold was 
shown bv lKamp fe van Zadelhofj l|200l|) . who calculated 
the gas temperature in tenuous disks around Vega-like 
stars. These disks have very low disk masses, of order a 
few Mq, and are optically thin to ultraviolet (UV) radia- 
tion throughout. It was shown that their gas temperature 
is generally very different from the dust temperature, and 
that this higher temperature significantly affects the inten- 
sity of gaseous emission lines, in t he most extreme c ases 
by an order of magnitude or more <|Kamp et alJ l2003). 

In this work, the gas temperature in the more massive 
disks (up to 0.1 Msun) around T-Tauri stars is examined, 
which are optically thick to UV radiation. Special atten- 
tion is given to the effects of dust settling and the influence 
of explicit gas temperature calculations on the properties 
of observable emission lines. 

The model used in this paper is described in Section 2. 
The resulting temperature, chemistry and emission lines 
are shown in Section 3. In Section 4 the limitations of our 
calculations are discussed, and in Section 5 our conclusions 
are presented. The details of the heating and cooling rates 
are in Appendices A and B, respectively. 

2. Calculating the gas temperature 

2.1. The PDR model 

Because of the symmetry inherent in the disk shape, it is 
usually assumed that disks have cylindrical symmetry, re- 
sulting in a 2-dimensional (2-D) structure. Because of the 
computational problems in solving the chemistry (partic- 
ularly the H2 and CO self shielding) and cooling rates in 
a 2-D formalism, however, these are calculated by divid- 
ing the 2-D structure into a series of 1-D structures (see 
Figure nj. The disk is divided into 15 annuli between 50 
and 400 AU. These 1-D structures resemble the photon- 



dominated regions (PDRs) f ound at the edges of molecu - 
lar clouds (for a review, see iHollenbach fc Tielenslll997l) . 
A full 2-D formalism is used to calculate the dominant 
heating rates due to the photoelectric effect on Polycyclic 
Aromatic Hydrocarbons (PAHs) and small grains. 

The 1-D PD R mod el described by 

iBlack fc van Dishoeckl lll987t) van Dishoeck fc Black 
(19881) and Ijansen et al .1 lll995h is used in this work. 
This code consists of two parts: the first part calculates 
the photodissociation and excitation of H2 in full detail, 
taking all relevant H2 levels and lines into account. 
It includes a small chemical network containing the 
reactions relevant to the formation and destruction of 
H2 to compute the H/H2 transition. The main output 
from this code is the H2 photodissociation rate as a 
function of depth and the fraction of molecular 
hydrogen in vibrationally-excited states. The second 
part of the program includes a detailed treatment of 
the CO photodissociation process, a larger chemical 
network to determine the abundances of all species, and 
an explicit calculation of all heating and cooling processes 
to determine the gas temperature. The latter calculation 
is done iteratively with the chemistry, since the cooling 
rates depend directly on the abundances of O, C + , C and 
CO. Each PDR consists of up to 130 vertical depth steps, 
with a variable step size adjusted to finely sample the 
important H— >H2 and C + ^C^CO transitions. 

To simulate a circumstellar disk, some modifications 
had to be made in the PDR code. Because of the higher 
densities, thermal coupling between gas and dust had to 
be included. While this is not an important factor in most 
PDRs associated with molecular clouds, it dominates the 
thermal balance in the dense regions near the midplane 
of a disk. A variable density also had to be included to 
account for the increasing density towards the midplane 
of the disk. As input to the code, the den sity structure as 
well a s the dust temperature computed bv lD'Alessio et all 
l|l999|) are used. In this model, the structure is calculated 
assuming a static disk in vertical hydrostatic equilibrium. 
The disk is considered to be geometrically thin, so ra- 
dial energy transport can be neglected. The model fur- 
ther assumes a constant mass accretion rate throughout 
the disk (a value of M = 10~ 8 Mq yr _1 is used), and tur- 
bulent viscosity given by the a prescription (a — 0.01). 
In the outer disk studied here, the contribution of accre- 
tion to the heating is negligible. The central star is as- 
sumed to have a mass of 0.5 Mq, a radius of 2Rq and a 
temperature of 4000 K. The disk mass is 0.07 Mq and 
ex tends to ~400 AU. Thi s same model has been u sed 
bv lAikawa et alJ l|2002j) and lvan Zadelhoff et all l)2003|) to 
study the chemistry of disks such as that toward LkCa 15 
and TW Hya. As shown in those studies, the results do 
not depend strongly on the precise model parameters. 

Further input to the PDR code is the strength of the 
UV field. The U V field is ba s ed on the spectrum for TW 
Hya obtained bvlCosta et ail l)2000|) (stellar spectrum B in 
Ivan Zadelhoff et all2003h : a 4000 K black body spectrum, 
plus a free-free and free-bound contribution at 3 x 10 4 K 
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and a 7900 K contribution covering 5% of the central star. 
The intensities of the resulting UV radiation incident on 
the disk surface at each annulus (with both a stellar and 
an interstellar comp onent) are calculated using the 2-D 
Monte Carlo code of Ivan Zadelhoff et alJ l|2003|) . The re- 
sulting UV intensities are then converted into a scaling 
factor Jtjv with resp ect to the stan dard interstellar radia- 
tion field as given bv iDrairiel (^978). These scaling factors 
Jtjv form the input to calculate the chemistry for each of 
the 15 vertical PDRs and range from ~ 1000 at the first 
slice at 6 3 AU to ~ 100 at the last slice at 373 AU. As 
shown by van Zadelhoff et aP (2003), the difference in the 
chemistry calculated with spectrum B and a scaled Drainc 
radiation field (spectrum A) is negligible. It is important 
to note that both radiation fields have ample photons in 
the 912-1100 A regime which are capable of photodisso- 
ciating H2 and CO and photoionizing C. 

The UV radiation also controls the heating of the gas 
through the photoelectric effect on grains and PAHs. The 
rates of these two processes are calculated using the 2-D 
radiation field at each radius and h eight in the disk com- 
puted bv lvan Zadelhoff et 2] §003) including absorption 
and scattering. 



2.2. Chemistry 



To ca lculate the chemistry, the network of Ijansen et alJ 
( 1995) is used. It contains 215 species consisting of 26 ele- 
ments (including isotopes) with 1549 reactions between 
them. The adopted abundances of the most important 
elements are shown in Table ^ the carbo n and oxygen 
abund ances are similar to those used by lAikawa et alJ 
ll2002h . The chemistry has been checked against updated 
UMIST databases, but only minor differences have been 
found for the species observed in PDRs. Since the temper- 
ature depends primarily on the abundances of the main 
coolants, the precise chemistry does not matter as long as 
the C + ^C^CO transition is well described. The adopted 
cosmic ray ionization rate is £ = 5 x 10~ 17 s _1 . 

Since the chemical timescales in the PDR layer are 
short (at most a few thousand yr) compared to the life- 
time of the disk, the chemistry is solved in steady-state. 
Only gas-phase reactions are considered; no freeze-out 
onto grains is included. Accordingly, the PDR calculations 
are stopped when the dust temperature falls below 20 K. 
Since this is also the regime where the gas and dust are 
thermally coupled (see §3.1), no explicit calculation of the 
gas temperature is needed. 

It should be noted that only the vertical attenuation 
is considered here in the treatment of the p hotorates used 
in the chemistry, in contrast to the work bv lAikawa et all 
(2002) who included attenuation along the line of sight 
to the star. This means that the photorates are gener- 
ally larger deeper into the disk. This is illustrated in 
Figure [2 where the dissociation rate of C2H is shown 
as functions of height. Comp arison with the rate found 
by Ivan Zadelhoff et all l(2003h shows that the dissociat- 
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Fig. 2. The photodissociation rate of C2H at a radius of 
105 AU. The solid line giv es the rate found in this w ork, 
the dotted line the rate bv lvan Zadelhoff et alJ l)2003j) . 

Table 1. Adopted gas-phase elemental abundances with 
respect to hydrogen. 



Element 


abundance 


D 


1.5 x 10" b 


He 


7.5 x 1(T 2 


c 


7.9 x 1(T 5 


N 


2.6 x 10~ 5 


O 


1.8 x 1(T 4 


s 


1.7 x 10~ 6 


PAH 


1.0 x 10~ 7 



ing radiation penetrates deeper (to a height of ~ 40 AU 
opposed to ~ 60 AU) into the disk. Given the over- 
all uncertainties in the radiation field and disk struc- 
ture these effects are not very s ignificant. Another d iffer- 
ence of our model with th ose of lAikawa et al.l ((2002) and 
Ivan Zadelhoff et all l|2003l) is the full treatment of CO pho - 
todissociation as given by Ivan Dishoeck fc Black! Jl988), 
including self shielding and the mutual shielding by H2. 

2.3. Thermal balance 

The gas temperature in the disk is calculated by solving 
the balance between the total heating rate T and the to- 
tal cooling rate A. The equili brium temperatur e is deter- 
mined using Brent's method ijPress et al.lll989h . Heating 
rates included in the code are due to the following pro- 
cesses: photoelectric effect on PAHs and large grains, cos- 
mic ray ionization of H and H2, C photoionization, H2 
formation, H2 dissociation, H2 pumping by UV photons 
followed by collisional de-excitation, pumping of [O I] 
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by infrared photons followed by collisional de-excitation, 
exothermic chemical reactions, and collisions with dust 
grains when Td us t > T gas . The gas is cooled by line emis- 
sion of C + , C, O and CO, and by collisions with dust 
grains when T dust < T gas . A detailed description of each 
of these processes is given in Appendix A. The thermal 
structure found in our PDR models has been checked ex- 
tensively against that computed in other PDR codes. 

There are two important assumptions in our calcula- 
tion of the thermal balance which may not necessarily be 
valid for disks compared with the traditional PDRs as- 
sociated with molecular clouds. The first is the use of a 
1-D escape probability method for the cooling lines. This 
limitation is further discussed in §4.1. The second is the 
assumption that the photoelectric heating efficiency for 
int erstellar grains also holds for g rains in disks. As shown 
bv lKamp fc van Zadelhofl lj200lF . this efficiency is greatly 
reduced if the grains have grown from the typical inter- 
stellar size of 0.1 pm to sizes of a few /im. While there 
is good observational evidence for grain growth for older, 
tenuous disks, the younger disks studied here are usually 
assumed to have a size distribution which includes the 
smaller grains. 



2.4. Dust settling 

In tur bulent disks, the gas and du st are generally well- 
mixed. IStepinski fe Valagead {l996) have shown that the 
velocities of dust grains with sizes < 0.1 cm are coupled 
to the gas. When the disk becomes quiescent, models pre- 
dict that dust particles are no longer supported by the gas 
and will start to settle towards the midplane of the disk 
l) Weidenschillind li 997j) . During the settling process, dust 
particles will sweep up and coagulate other dust particles, 
thus leading to grain growth. Since larger particles have 
much shorter settling times than small particles, coagula- 
tion accelerates the settling process. 

Observational evidence for dust growth and settling 
in disks includes the near-IR morphology of the edge- 
on disks. For exam ple, for the 114-526 disk in Orion, 
iThroop et ail (|200l[) and lShuping et all l|2003j) show that 
the large dust grains in this disk are concentrated in the 
midplane. The SEDs of some o f the T-Tauri stars ex- 
amined bv lMivake fc Nakagaw also indicate that 
dust settling takes place, as well as the stat i stics o f obser- 
vations of edge-on disks by iD'Alessio et alJ l|l999|l . 

In our model, dust settling is simulated by varying the 
gas/dust mass ratio. In this paper, both models with a 
constant m gas /rodust ( "well- mixed" ) and a variable value 
of 

^gas /^dust ("settled") are presented. In the well-mixed 
case, the gas/dust mass ratio is taken to be 100 through- 
out. For the settled model, the value of 100 is kept near 
the midplane of the disk, while a value of 10 4 is used in 
the surface regions. The precise value adopted at the sur- 
face is not important, as long as it is high enough that 
photoelectric heating is no longer significant in this layer. 
The ratio is varied linearly with depth in a narrow tran- 
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Fig. 3. Settling times in years for 0.1 /Ltm-sized dust parti- 
cles. The dashed lines give the boundaries between which 
the rrigas/mdust ratio is interpolated from a value of 10 
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sition region. The boundaries of the transition region are 
defined as z% ±/i/10, where z 6 is the height where the set- 
tling time of dust particles is 10 6 years (this method thus 
simulates a disk that is 10 6 years old) and h is the height 
of the DAlessio disk, defined as the height where the gas 
pressure is equal to 7.2 x 10 5 K cm -3 . 

The settling time t S ettie is determined as follows: 

7. 

^settle = 



^settle 



where z is the height above the midplane. The settling 
velocity ifeettie is determined by the balance between the 
vertical component of the gravitational force on the dust 
grains and the drag force the grain experiences as it moves 
through the gas. For the drag for ce the subsonic limit of 
the formulae by IBerruverl ljl99l|) is used, resulting in a 
settling velocity of 



^settle 



y/n G M zap 
2 r 3 p n H «th. 



crmal 



with G the gravitational constant, M the stellar mass, a 
the grain radius (a typical value of 0.1 pva is used), p the 
density of the grain material (2.7 g cm -3 ), r = y/ z 2 + R 2 
the distance from the star, p the mean molecular weight 
(2.3 proton masses) and ^thermal the thermal velocity of 
the gas particles (v thermal = \/2 k T/p). The settling times 
found by this method are displayed in Figure |3 

It is further assumed that the smallest particles rep- 
resented by PAHs remain well-mixed with the gas in all 
models considered here. Observations show evidence for 
PA H emission from dis ks, at least around Herbig Ae stars 
fe.g. lMeeus et al.ll200l|) . whereas detailed models indicate 
that their settling times are muc h longer than the ages of 
the disks l| Weidenschillind 1 1 997h . This means that PAH 
heating is still at full strength in the upper layers of the 
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Fig. 4. The temperature structure of the disk. The colorscale and the white contours give the gas temperature, the 
black contours give the dust temperature for the well-mixed model (a) and the dust settling model (b). The dominant 
cooling rates are also shown for the well-mixed (c) and settled (d) models, as well as the dominant heating rates in 
(e) and (f). 



disk while photoelectric heating on large grains is sup- UV radiation with wavelen gths shorter than 1500 A which 
pressed in the case of dust settling. The PAHs are also dissociates molecules (e.g. iLi fc Greenberdll99^) . Thus, 
responsible for approximately half of the absorption of much of the UV radiation which affects the chemistry is 
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Figures (c) and (d) give the cooling rates at this radius, where the solid line denotes cooling by gas-dust collisions, the 
dashed line [C i] cooling, the dotted line CO cooling, the dash-dot line [C il] cooling and the dash-triple dot line [O i] 
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rate due to gas-grain collisions. 



Jonkheid et al.: The gas temperature in flaring disks 



7 



still absorbed in the upper layers even when large dust 
particles are settling. In our models, this is taken into ac- 
count by adopting an effective Ay in the calculation of 
the depth-dependent photodissociation rates. Specifically, 
the gas/dust parameter x g( j is defined as the actual value 
of the ra gas /rad U st divided by the interstellar value. Thus, 
Xgd varies between 1 and 100. The visual extinction then 
becomes 

A = ( Nr \ x — 
V ^1.59 10 21 cm- 2 J X x gd 

The effective extinction used for photorates at UV wave- 
lengths 

Av ' cS = (l.591CFcm-0 X + 

Thus the dissociating radiation is still reduced (albeit at 
half strength) in areas where large grains have disappeared 
completely. 

The scaling of the various heating and cooling rates in 
the settling model is described in Appendix A. In particu- 
lar, it should be noted that although the H2 grain surface 
formation rate is reduced, H2 formation through the reac- 
tion PAH : H + H — * PAH + H 2 is included in the chemi- 
cal network. Since this reaction has a rate comparable to 
the H2 formation rate on large grains the abundance of 
H2 is lowered by only a factor of 2 in the upper layers 
if the settled disk. If H2 formation through PAHs is ex- 
cluded, the H — > H 2 transition occurs deeper in the disk, 
but the effect on the C + — ► C — > CO transition is small. 
The overall gas density structure is kept the same as in 
the well-mixed model, as is the dust temperature distri- 
bution. The latter assumption is certainly not valid, but 
since gas-dust coupling only plays a role in the lower lay- 
ers, this does not affect our results. 

3. Results 

3.1. Temperature 

The results of the temperature calculations are presented 
in Figure^ for the entire disk, and in Figure [S] for a ver- 
tical slice through the disk at 105 AU. It can be seen 
that the vertical temperature distri bution resembles that 
of a " normal" PDR (see for examnle lTielens fc Hollenbachl 
1985): the gas temperature is much higher than the dust 
temperature at the disk surface, and both temperatures 
decrease with increasing visual extinction. In annuli close 
to the central star the surface temperatures go up to 1000 
K. The precise temperature in these annuli is uncertain 
and depends on the adopted molecular parameters. For 
example, when the H* vibrational de-excitat ion rate coeffi- 
cients given bv lTielens fc Hollenbachl l|l985|) are used, the 
rate increases steeply with increasing temperature with 
the cooling rate unable to keep up, resulting in a numeri- 
cal instabjhty_^ntlie code. W hen the rate coefficients given 
bv lLe Bourlot et alJ l)l999l) are used, a stable temperature 
can be found. 



10° 






10- 1 




-I 


* 10-2 




_ 








\ 10- 3 
\ 1 w 












10- 4 






10- 5 






If)" 6 







200 400 600 800 1000 1200 
T (K) 

Fig. 6. The normalized distribution of the temperature 
over total disk mass. The solid line shows the dust tem- 
perature, the dashed line the gas temperature of the well- 
mixed model, and the dotted line the gas temperature of 
the settled model. 

The main difference between the temperature struc- 
tures in disks and general molecular clouds is caused by 
the increasing density in disks: in standard PDRs the den- 
sity generally stays uniform, causing the gas temperature 
to fall below the dust temperature at high optical depths. 
In disks, the increasing density causes the coupling be- 
tween gas and dust to dominate the thermal balance in 
the deeper layers, so that the gas temperature becomes 
equal to the dust temperature. 

The normalized distribution of the gas temperature 
over the disk mass between 63 AU and 373 AU of the 
central star is shown in Figure El The calculations are 
performed down to a dust temperature of ~ 20 K. It can 
be seen that even though the bulk of the gas mass is at 
the same temperature as the dust, a considerable fraction 
(~ 10%) in the surface layers is at higher temperatures, 
especially in the case of dust settling. Since most of the 
molecular emission arises from this layer, the difference 
is significant. It can also be seen that the settled model 
contains the largest amount of hot gas, even though the 
maximum temperature is higher in the well-mixed case. 

The above figures illustrate that T gas = Tdust is in- 
valid in the upper layers of the disk. This can also be seen 
in the individual heating and cooling rates in Figure 
gas-dust collisions dominate the heating/cooling balance 
only deep within the disk. In the upper layers the heating 
due to the photoelectric effect on large grains and PAHs 
is so large that cooling of the gas through collisions with 
dust grains is not effective anymore, and atomic oxygen 
becomes the dominant cooling agent. This shifts the equi- 
librium to higher gas temperatures. 

In the case of dust settling, the photoelectric effect 
on large grains and the H2 formation heating rates are 
suppressed by a factor of 100. The PAH heating is still 
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Fig. 7. Chemical abundances for the well-mixed (left column) and settled (right column) models. In figures (a) and 
(b), the density profile of the disk is shown, as well as the Ay=l surface denoted by the dashed line. In figures (c) and 
(d), the surface of the disk is shown as a dotted line, the H/H2 transition as a dashed line, and the O/O2 transition 
as a solid line. The shaded area denotes the C+/C/CO transition. 



effective, resulting in a temperature that is ~ 50% lower 
than in the well-mixed model. In this model there is also 
less absorption of UV radiation in the upper layers, result- 
ing in higher temperatures deeper in the disk compared 
to the well-mixed model. When the abundance of large 
grains rises, all this UV radiation is available for photo- 
electric heating, resulting in a peak in the gas tempera- 
ture. This extra UV is rapidly absorbed here, causing a 
sudden drop in the temperature as gas-dust collisions be- 
come the dominant process in the thermal balance. If the 
PAHs are removed from the model, thus simulating a disk 
where these small grains have disappeared, the tempera- 
ture in the outer regions drops even further, but the overall 
shape of the temperature structure does not change much. 

3.2. Chemistry 

In Figures and |H1 the chemical structure of the disk is 
presented. The chemistry follows that of an ordinary PDR: 



a region consisting mainly of atomic H in the outer lay- 
ers, with a sharp transition to the deeper molecular layers. 
Self-shielding of H2 quickly reduces the photodissociation 
rate, so the deeper layers consist mainly of molecular hy- 
drogen. The principle carbon-bearing species follow a simi- 
lar trend, consisting mainly of C + in the upper layers, with 
a transition to C, and later CO, at larger optical depths. 
In the deeper layers CO is also protected from dissociation 
by self-shielding and by shielding by H2, as several disso- 
ciative wavelengths overlap for these molecules. Because 
CO is much less abundant than H2, absorption by dust 
also plays an important role in decreasing its dissociation 
rate. This causes the carbon in the disk to be mostly ionic 
in the upper layers, and mostly molecular (in the form 
of CO) deepe r in th e dis k. Compared with the mod els of 
lAikawa et all ((20021) and Ivan Zadelhoff et alJ l|2003h . our 
C — * CO transition occurs somewhat deeper into the disk, 
at z —40 AU rather than 60 AU for the i?=105 AU annu- 
lus. This is due to the geometry in the 1+1-D model, which 
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Fig. 8. The vertical distribution of several chemical species in the disk at 105 AU for the well-mixed (left column) 
and settled (right column) models. Figures (a) and (b) show the abundances of H and H2; figures (c) and (d) of C + , 
C and CO. 



allows dissociating radiation to penetrate more deeply into 
the disk, as well as our different treatment of CO photodis- 
sociation. 

In the case of dust settling the C + /C/CO transition 
occurs slightly deeper in the disk. This is explained by the 
higher photodissociation rates deeper in the disk due to 
the decreased absorption of UV radiation by large dust 
grains. The PAHs in the outer regions still absorb a sig- 
nificant fraction of the UV, so the effect is not very large. 
Since the H/H2 transition is determined by H2 self shield- 
ing and not dust absorption, the position of this transition 
does not change if the dust settles. 

The chemistry has also been calculated for a model 
where T gas = Tdust- It is found that the different temper- 
ature has little effect on the abundances of those species 
important in the thermal balance. The chemistry in the 
surface layers of disks is mostly driven by photo-reactions 
and ion- molecule reactions, both of which are largely in- 
dependent of temperature. 



3.3. Line intensities 

The main effect of the high gas temperatures is on the 
intensities and shapes of emission lines arising from the 
warm surface layers. It is expected that the intensities 
of lines tracing high temperatures (such as the [0 1] fine- 
structure lines, the H2 rotational lines, and the higher CO 
rotational lines) will increase due to the higher tempera- 
tures, and the larger overall amount of warm gas. 

The model line p rofiles are created using the 2-D 
Monte Carlo code by iHogerheiide fc van der Takl (2000) 
assuming a Keplerian velocity field and an inclination of 
60°. The results for the [Cl], [Oi] and [C11] fine-structure 
lines and the rotational lines of CO are shown in Figure^ 
It can be seen that the high gas temperatures have a sig- 
nificant impact on the intensities, and particularly on the 
[Oi] and [C11] fine-structure lines. The difference in the 
locations in the disk where these lines are predominantly 
excited is reflected in the different shapes of the lines. The 
lines tracing hot gas all display a double-peaked structure, 
due to their excitation in the innermost parts of the disk. 
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Fig. 9. The emission lines of C, O, C + and CO for three scenarios: T gas = Tdust (dotted line), and T gas calculated 
explicitly in the well-mixed (solid line) and settled (dashed line) model. Units are in Kelvins, assuming that the disk 
fills the beam precisely; this corresponds to a beam size of 5" when the disk is at 150 pc. 



It can also be seen that the higher gas temperature has 
little effect on the lower rotational lines of CO and the 
lowest [C i] fine-structure line. These lines trace the cold 
gas, which is present in the deeper layers of the disk in all 
models. 



4. Discussion 

4.1. Limitations of the 1+1-D model 

While circumstellar disks are inherently 2-dimensional 
structures, a complete 2-D treatment of the radiative 
transfer, chemistry and thermal balance is at present too 
time-consuming to be practical. The 1+1-D model circum- 
vents this problem by drastically simplifying the geometry 
to a series of 1-D structures. The main disadvantage of the 
1+1-D model is its poor treatment of radiative transfer: 
in principle, only transfer in the vertical direction is in- 
cluded. In a pure 1+1-D geometry the stellar radiation 
is forced to change direction when it hits the disk's sur- 
face and thus affects the chemistry calculation. This prob- 
lem does not apply to the photoeletric heating rates of 
PAHs and large grains since they were calculated using 



a UV field calculated with the 2-D Monte Carlo code by 
Ivan Zadelhoff et all l)2003(l . 

The 1+1-D geometry also affects the escape probabil- 
ities of the radiative cooling lines, which have to travel 
vertically instead of escaping via the shortest route to the 
disk's surface. Thus the escape probabilities calculated by 
the 1+1-D model are generally too low. 

The overall effect of the 1+1-D geometry on the disk 
temperature is subtle: the very upper layers of the disk 
are optically thin regardless of the adopted geometry, so 
the temperature and chemistry results are expected to be 
correct there. In the dense regions near the midplane, the 
dust opacities are so large that there will be little effect of 
the geometry either. Also, due to the strong thermal cou- 
pling between gas and dust, the gas temperature is equal 
to the dust temperature there. The greatest uncertainties 
are in the chemistry, in the intermediate regions between 
the surface and the midplane. Since the thermal balance 
is tied to the chemistry, it is expected that the largest 
uncertainties in the gas temperature occur also in this re- 
gion. It should be noted that many secondary effects play 
a role here, including the precise formulation of the cool- 
ing rates, treatment of chemistry, self-shielding, etc. The 
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Fig. 10. Vertical distribution of gas and dust tempera- 
tures in the disk at a radius of 105 AU. The solid line 
gives the gas temperature when gas-dust collisions are ig- 
nored, the dashed line gives the gas temperature when 
these collisions are included in the thermal balance. The 
dotted line gives the dust temperature. The inset shows 
the region around T gas = 20 K. 

current models should be adequate to capture the main 
characteristics and magnitude of all of these effects. 

4.2. Effects of gas-dust coupling 

It can be seen in Figure HUI that thermal coupling between 
gas and dust has a small effect on the gas temperature 
in the upper layers of the disk: the gas temperature is so 
high and density comparitively low that gas-dust collisions 
are an ineffective cooling mechanism. Deeper in the disk 
the coupling becomes stronger and eventually comes to 
dominate the thermal balance. The overall effect on the 
resulting temperature is small because gas and dust al- 
ready have similar temperatures even when the coupling 
is ignored. 

4.3. Effects of the radiation field 

The radiation field used in th is work is ass umed to have 
the same spectral shape as the lDrainel ((1278) field with in- 
tegrated inte nsity matched to the observed radiation field 
of TW Hya l|Costa et alJl2000jh iBergin et all l|2003h have 
shown that this treatment may overestimate the amount 
of radiation in the 912-1100 A regime where H 2 and CO 
are dissociated, since a significant fraction of the UV flux 
is in emission lines (particularly Ly a). Thus the dissoci- 
ation rates of H2 and CO, as well as the C ionization rate 
may be too high in this work. Other molecules such as 
H2O and HCN, however, have significant cross sections at 
hya. 

If the radiation field of IBergin et al.l l|200.ll) is used, 
the chemistry of the relevant cooling species will become 



more similar to that found using s pectru m C (a 4000 K 
blackbody) in Ivan Zadelhoff etafl |£o03). Because of the 
different chemistry, it is also expected that the temper- 
ature profile will change: CO will become the dominant 
carbon-bearing species at larger heights, where it can cool 
the gas more efficiently. It is therefore expected that the 
temperature in the intermediate layers will drop with re- 
spect to the results presented here. The heating rates are 
not sensitive to changes in the chemistry, so they will not 
change much. 



5. Conclusions 

The main conclusions of our calculations are: 

— The gas temperature is much higher than the dust tem- 
perature in the optically thin part of the disk, while 
the two temperatures are the same in the part of the 
disk which is optically thick to UV radiation and where 
most of the disk's material resides. 

— The temperatures in the optically thin part have a no- 
ticeable effect on the excitation of higher-lying emis- 
sion lines. The higher temperatures should also affect 
the disk's structure: since the structure is calculated 
assuming hydrostatic equilibrium, the high temper- 
atures in the surface regions will blow-up the disk, 
which will achieve hydrostatic equilibrium at a greater 
height. This in turn will lead to a more efficient out- 
gassing of the disk. An initial d i scussi on of this effect 
is given in lKamp fc Dullemondl ((2004). 

— The chemistry in the disk does not change much when 
the gas temperature is increased in the explicit calcu- 
lation. 

— Dust settling is found to have a great impact on the gas 
temperature in the disk. The gas temperature in disks 
where settling is taking place is found to be lower, and 
decreases less steeply than in well mixed disks. Due to 
the lower temperatures, the intensities of higher-lying 
lines are lower than in well mixed disks. 

— Dust settling does not greatly affect the disk's chem- 
istry. This is because PAHs are assumed to stay well- 
mixed with the gas when large grains are settling. As 
a result much of UV radiation important for chemistry 
is still absorbed in the surface layers of the disk in the 
settled model. The C+/C/CO and the 0/0 2 transi- 
tions thus occur only slightly lower in the disk. 

— Dust settling increases the intensity of atomic and 
molecular lines, and has a subtle influence on the line 
shapes due to the smaller amount of hot gas, especially 
at small radii. 
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Appendix A: Heating processes 

In the following, the various processes included in the cal- 
culation of the total heating rate are discussed. Also, their 
scaling with m gas /77id U st is indicated. 

Photoelectric heating: the energetic electrons released 
by absorption of UV photons by dust grains contribute 
significantly to the total heating rate in the surface of the 
disk. For the h eating rate due to the photo electric effect 
the formula bv lTielens fc Hollenbachl l|l985h is used: 

r PE = 2.7 x 10 - 25 5 UV 5 d n H Y I vv 
x [(1 — x) 2 /x + Xk(x 2 — l)/x 2 ] ergcm~ 3 s _1 
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with S uv = 2.2, 6 d = 1.25 and Y = 0.1. The strength of 
the radiation f ield Ittv is calculated with the 2-D Monte 
Carlo code bv Ivan Zadelhoff et alJ l)2003|) . x is the grain 
charge parameter, found by solving the equation 



+ (xk - Xd + l)x 2 - 7 







where x = Eq/Eh, Xk — kT/En, Xd — Ed/ Eh and 
7 = 2.9 x 10~ 4 YV T ri\jyn~ 1 ; Eh is the ionization po- 
tential of hydrogen, Ed is the ionization potential of a 
neutral dust grain (6 eV is assumed here) and Eq is the 
grain potential including grain charge effects. This rate is 
scaled with m^ust /mgas to account for the varying gas to 
dust ratio in models where dust settling was included. 

PAH heating: analogous to photoelectric heating, the 
photoelectrons from PAH ionization also contribute to 
the heating of the gas . The heating rate described in 
iBakes fc Tielensl i|l994|) is used, using PAHs containing 
Nq carbon atoms (30 < < 500). The radiation field 
used to obtain this heat ing rate is calculated with the 
2-D Monte Carlo code bv lvan Zadelhoff et alJ l|2003|) . 

C photoionization: it is assumed that each electron re- 
leased by the photoionization of neutral C delivers approx- 
imately 1 eV to the gas. The heating rate then becomes 

T c ion = 1.7 x 10~ 12 Rc io„ n(C) erg cm" 3 s" 1 

where Rc ion is the C ionization rate. 

Cosmic rays: another source of energetic electrons is 
the ionization of H and H2 by cosmic rays. It is assumed 
that about 8 eV of the electron's energy is used to heat 
the gas in the case of H 2 , and 3.5 eV in the case of H. 
The H2 contribution is scaled to include the ionization of 
He. A cosmic ray ionization rate of £cr = 5 x 10" 17 s _1 is 
used. The heating rate then becomes 



Tcr = [2.5 x 10~ n n(H 2 ) + 5.5 x 10" 
x Ccr erg cm -3 s -1 



.(H)] 



H2 formation: of the 4.48 eV liberated by formation of 
H2 on dust grains, it is assumed that 1.5 eV is returned 
to the gas. The corresponding heating rate is: 

Th 2 form = 2.4 X 10" 12 i?H 2 form ™H erg Cm" 3 S^ 1 

where i? H2 form = 3 x 10~ 18 VTnjj n(H)/(l + T/T stick ) 
is the H2 formation rate, with T st i c k=400 K. This rate is 
scaled with TOdust/TO gas in models where dust settling is 
included. 

H2 dissociation: when H2 is excited into the Lyman 
and Werner bands, there is about 10% chance that it will 
decay into the vibrational continuum, thus dissociating 
the molecule. Each of the H atoms created this way carries 
approximately 0.4 eV. This gives: 



where i?n 2 d i aB ^ s * ne H2 photodissociation rate. 

H2 pumping: electronically excited H2 that does not 
disscociate decays into bound vibrationally excited states 
of the electronic ground state. It is assumed that 2.6 eV is 
returned to the gas by collisional de-excitations, resulting 
in: 

r H2 pump - 4.2 x HT 12 [n(H) 7H + n(H 2 ) 7 H 2 ] 
x nCH.^) ergcm~ 3 s 

where 7h and 7h 2 are the deexcitation rate coefficients 
through collisio ns by H and H?, r espec tively. These rates 
are taken from lLe Bourlot et al.1 l)l999j) . H?; is calculated 
explicitly in our models in the UV excitation of H2. 

Chemical heating: in principle every exothermic reac- 
tion contributes to the heating of the gas. The reactions 
considered here (and the energies released) are the disso- 
ciative recombination of H^ (in which 9.23 eV is realeased 
by dissociating to H+H 2 and 4.76 eV by dissociating to 
H+H+H), HCO+ (7.51 eV) and H 3 0+ (6.27 eV), and the 
destruction by Hc+ ions of H 2 (6.51 eV) and CO (2.22 
eV). 

Gas-dust collisions: collisions between gas and dust 
will heat the gas when the dust temperature is higher 
than the gas temperature, and will act as a cooling term 
when the dust tempe rature is lower than the gas t empera- 



ture. The formula bv lBurke fc Hollenbachl l|1983D is used, 
assuming an average of (<7 gr n gr ) = LSnncm" 1 and an 
accomodation coefficient olt = 0.3. The resulting rate is 
scaled with the local value of TOdust /w gas in models where 
dust settling was included: 

r gd = 1.8 x IO- 33 n| VT (T dust - T gas ) erg cm- 3 s" 1 



Appendix B: Cooling processes 

The gas is cooled through the fine-structure lines of C + , 
C and O, and the rotaional lines of CO. For each species 
the level population was calculated assuming statistical 
equilibrium: 

dri{ ^ — ^ — ^ 
— = 2^ rij Rj^i - Ri^ 







where rii is the density of a given species in energy state 
i, and Ri-,j is the total rate of level i to level j. For all 
species, collisional excitations and de-excitations are taken 
into account, as well as spontaneous emission and absorp- 
tion and stimulated emission through interaction with the 
cosmic microwave background (CMB) and the dust in- 
frared background. For C + and O, UV pumping of the 
fine-structure lines is also included. Once the populations 
ar e known, the cooling; rate can be found using the formula 
bv lTielens k Hollenbachl l)l985|) : 



Fh, 



= 6.4x10 i3 7i(H 2 )i?H 2 diss erg cm 3 s 



S{vjj) -P(vjj) 



14 



Jonkheid et al.: The gas temperature in flaring disks 



In this formula, rii is the density of species X in energy 
state i, Aij is the Einstein A coefficient for the transition 
i —y j, Vij is the corresponding frequency, and P C sc( T ij) is 
the escape probability at optical depth ry. S(vij) is the 
local source function: 

2K" (9i*i 



%) = -rP-i 

where gi is the statistical weight of level i. P(vij) is the 
intensity of the background radiation: 

P(fij) = B Vij (T C mb) + T dnst B Vi] (T ) 

where Xcmb is the temperature of the cosmic background 
(2.726 K) and To the characteristic temperature of the 
dust at the surface of the disk. In the calculation of the 
collision rates of all cooling species considered here, the 
ortho/para ratio for H2 was assumed to be in thermal 
equilibrium with the local gas temperature. 

cooling: O contributes to the cooling of the gas in the 
surface layers of the disk through its fine-structure lines 
at 63.2 /jm and 145.6 /im. The Einstein A coefficients for 
these lines are 8.87 x 10~ 5 s _1 and 1 .77 x 10~ 5 s -1 , respec- 
tively . C ollisions with electrons es fe Nussbaumerl 



1984 - H l|Launav fc Rouefi1ll9774 and H9 l|Jaauet et al 



19921 were included. O can also act as a heating agent if 
infrared pumping is followed by collisional de-excitation. 
If that happens the cooling rate becomes negative and is 
treated as a heating rate by the code. This phenomenon 
did not occur in the calculations presented in this paper. 

C + cooling: C + contributes to the cooling of the gas 
through its fine-structure line at 158 /im. An Einstein A 
coeffic ient of 2.29 x 10~ 6 s~ 1 is use d. C ollisions with elec- 
trons jHaves fc Nussbaumerl 1 19841) H llLaunav fc RouefJ 
Il977bh and H9 ijFlower fc Launavlll977j) are included in 
calculating the level populations. 

C cooling: C cools the gas through its fine-structure 
lines at 370 /im and 609 /im. The Einstein A coefficients 
used are 2.68 x 10~ 7 s ~ 1 and 7.93 x 10~ 8 s~ 1 ] respec- 
tivelv. Collisions with H ( Launav fc Rouefjll977a|) and H 2 



ijSchroder et aljfl991^ are included. 

CO cooling: CO is the main coolant of the dense, 
shielded regions deep in the disk through its rota- 
tional lines. The 20 lowest rotational transitions of 
CO are included in the model. The Einstein A coeffi- 
cien ts used in this work ar e ident ical to those adopted 
bv iKam p fc van" Zadelhofj ll200lh Collisions with H 
l|Warin et al.lll996|) and H9 l)Flowerll20oH) are taken into 
account to calculate the populations. 



